PROGRAM GIS2BDY

  IMPLICIT NONE

  REAL           :: alatb, alonl, alatt, alonr
  REAL           :: dlat,  dlon,  plat,  plon
  INTEGER        :: nlat,  nlon,  kret,  ii, jj
  INTEGER        :: znot0, znot1, hgts0, hgts1, luse0, luse1

  CHARACTER(80)  :: fname

  INTEGER(4),    ALLOCATABLE :: LUSE(:,:)    ! land-use 11 categories
  INTEGER(4),    ALLOCATABLE :: ZNOT(:,:)    ! roughness length (cm)
  INTEGER(4),    ALLOCATABLE :: HGTS(:,:)    ! terrain height (m)

!-------------------------------------------------------------------
! write BDYFILES identification file ASCDATA.CFG, file also identifies
! default values for land-use and roughness length when the domain of 
! the data defined does not cover the extent of the meterological data

  OPEN(10,FILE='ASCDATA.CFG',ACTION='WRITE')

! lower left corner of the file
  ALATB= 33.0
  ALONL=135.0
  WRITE(*,*)'Lower left corner of file (lat,lon):',ALATB, ALONL
  READ (*,*) ALATB, ALONL

! lower left corner of the file
  ALATT= 42.0
  ALONR=143.0
  WRITE(*,*)'Upper right corner of file (lat,lon):',ALATT, ALONR
  READ (*,*) ALATT, ALONR

! grid spacing
  DLAT=0.1
  DLON=0.1
  WRITE(*,*)'Grid point spacing (dlat,dlon):',DLAT, DLON
  READ(*,*)  DLAT, DLON

! number of points in each direction
  NLAT=(ALATT-ALATB)/DLAT
  NLON=(ALONR-ALONL)/DLON
  WRITE(*,*)'Number of grid points (nlat,nlon)',NLAT,NLON
  IF(NLON.GT.1440)THEN
     WRITE(*,*)'Exceeding longitude dimension ... recompile'
     STOP
  END IF

! content default values
  LUSE1=7
  WRITE(*,*)'1-urban      2-farming      3-dry range 4-deciduous'
  WRITE(*,*)'5-coniferous 6-mixed forest 7-water     8-desert' 
  WRITE(*,*)'9-wetlands  10-mixed range  11-rocky'
  WRITE(*,*)'Land-use:',LUSE1
  READ(*,*)  LUSE1    

  ZNOT1=10
  WRITE(*,*)'Roughness length (cm):',ZNOT1
  READ(*,*)  ZNOT1

  HGTS1=0
  WRITE(*,*)'Terrain height (m):',HGTS1
  READ(*,*)  HGTS1

  WRITE(10,*)ALATB,ALONL
  WRITE(10,*)DLAT,DLON
  WRITE(10,*)NLAT,NLON
  WRITE(10,*)LUSE0
  WRITE(10,*)ZNOT0*0.01  ! in meters, although file content is cm
  WRITE(10,*)"'./'"
  CLOSE(10)

  ALLOCATE (LUSE(nlat,nlon), ZNOT(nlat,nlon), HGTS(nlat,nlon))

! initialize array to defaults
  LUSE=LUSE1
  ZNOT=ZNOT1
  HGTS=HGTS1

!------------------------------------------------------------------------
! open and read file with new land-use, roughness length, and terrain
! height values organized by record with each record identified by the
! center latitude-longitude location of the grid cell. This code section
! should be replaced user supplied content. Land-use is given in eleven
! categories, roughness length in centimeters, and terrain height in meters.

  WRITE(*,*)'Enter file name of new surface boudary data to be remapped:'
  READ(*,'(A)') fname
  OPEN(20,FILE=FNAME)

  kret=0
  DO WHILE (kret.EQ.0)
     READ(20,*,IOSTAT=kret) plat, plon, LUSE0, ZNOT0, HGTS0
     IF(kret.NE.0)THEN
        ii=1+NINT((plon-alonl)/dlon)
        jj=1+NINT((plat-alatb)/dlat)
    
        IF(ii.GE.1.AND.ii.LE.nlon.AND.jj.GE.1.AND.jj.LE.nlat)THEN
           LUSE(jj,ii)=LUSE0
           ZNOT(jj,ii)=ZNOT0
           HGTS(jj,ii)=HGTS0
        END IF
     END IF
  END DO
  CLOSE(20)
!-------------------------------------------------------------------------


! the new content is written to an ASCII output files in a simple 4-character
! integer format, each record, by latitude, starting from north, and with
! individual data values for each longitude point from west to east.

  OPEN(30,FILE='LANDUSE.ASC',ACTION='WRITE')
  DO JJ=nlat,1,-1
     WRITE(30,'(1440I4)') (LUSE(jj,ii),ii=1,nlon)
  END DO
  CLOSE(30)

  OPEN(40,FILE='ROUGLEN.ASC',ACTION='WRITE')
  DO JJ=nlat,1,-1
     WRITE(40,'(1440I4)') (ZNOT(jj,ii),ii=1,nlon)
  END DO
  CLOSE(40)

  OPEN(50,FILE='TERRAIN.ASC',ACTION='WRITE')
  DO JJ=nlat,1,-1
     WRITE(50,'(1440I4)') (HGTS(jj,ii),ii=1,nlon)
  END DO
  CLOSE(50)

END PROGRAM gis2bdy